#include <CRasterizer.h>


bool galaxySort (CGalaxy* a,CGalaxy* b) { 
  return ((a->params.position-a->params.camera).Length()> (b->params.position-b->params.camera).Length());
}


void CRasterizer::Prepare() {
    sort (galaxies.begin(), galaxies.end(), galaxySort); 

  }

#ifndef STRIP

void CRasterizer::Render() {
  AllocateData();

  //  mpi = new CMPI(MPI::COMM_WORLD.Get_size(), &components);

  mpi = new CMPI(MPI::COMM_WORLD.Get_size(), &image);
 
  RP.idx = 0;
  //  cout << "Rendering: [";
  
  // exit(1);
  int divide = RP.size*RP.size /(RP.no_procs);
  usleep(10);
  for (int i=1;i<RP.no_procs;i++) {
    RP.idx = i*divide;
    RP.idxCount = divide;
     mpi->SendToSlave(RP,i);
  }
  
  RP.idx = 0;
  RP.idxCount = divide;

  RenderPart(true);
  
  mpi->Test();
  //  cout << "]" << endl;

}

void CRasterizer::Master() {
    cout << "Rendering to size "<< RP.size << endl;
    Render();
    
    RP.idx = -1; // terminate

    for (int i=1;i<RP.no_procs;i++) {
      mpi->SendToSlave(RP,i);
    }

}


void CRasterizer::Slave() {
  CReturnValue v;
  double* data = NULL;

  while (true) {
    
    RP.dataBuffer = (double*) &RP;
    
    MPI::COMM_WORLD.Recv(RP.dataBuffer, RP.bufferSize, MPI::DOUBLE_PRECISION , 0, MPI::ANY_TAG);   

    if (RP.idx==-1) {
      return;
    }

    int w = RP.size;

    if (data==0)
      data = new double[(int)RP.idxCount];

    for (int i=0;i<RP.idxCount;i++)
      data[i] = 0.0;
    

    RenderPart(false);

    for (int i=0;i<RP.idxCount;i++) {
      data[i] = image.Data[i + (int)RP.idx].height.real();
    }

    MPI::COMM_WORLD.Send( data, RP.idxCount, MPI::DOUBLE_PRECISION , 0, 0);
    delete[] data;
    data = 0;
  }
  
}



#endif
void CRasterizer::AllocateData() {
  if (image.Data == 0) {
    image.Create(RP.size,RP.size);
  }
  image.Clear();
}



CGalaxy* CRasterizer::AddGalaxy(string file, CVector position, CVector orientation, double iscale) {

  if (spectra.size()==0) {
    throw new string("Error spectra needs to be loaded before galaxies are added");
  }

    CGalaxy* g = new CGalaxy();
    g->params.perlin = RP.perlin;
    g->params.orientation = orientation;
    g->params.position = position;
    g->params.intensityScale = iscale;
    g->loadSpectra(spectra);

    // randomize spectra
 
    g->InitializeFromFile(file, RP);

    if (g->components.size()==1) {
      // elliptical galaxy
      //g->params.bulgeAxis.y += (rand()%1000/1000.0 - 0.5)*1.0;
    }
    if (RP.randomize_spectra==1)
      g->RandomizeSpectra();

    galaxies.push_back(g);
    return g;
}


void CRasterizer::RenderPart(bool isMaster) {
  AllocateData();
  CCounter* cnt = 0;
  cnt = new CCounter(RP.idxCount*RP.no_procs, "Rendering image", false);

  for (int i=0;i<galaxies.size();i++) 
    galaxies[i]->setupComponentWavelength(RP.wavelength);



  for (int idx = RP.idx; idx<RP.idxCount + RP.idx;idx++) {
    
    cnt->Tick();
    setupCamera(idx);

    if (RP.lensing==1) 
      setupLensing(idx);

    renderPixel(idx);

  }
}

void CRasterizer::LoadRenderlist(string fname) {
  CUtil::verify_file(fname);
    
  fstream f(fname.c_str(), ios::in);
  if (!f)
    throw string("Could not load file "+fname);
  char s[2000];
  string er = "";
  while(!f.eof()) {
    f.getline(s,2000);
    vector<string> tok;
    CUtil::Tokenize(s, tok," ");
    if (tok.size()==8) {
      string name = "galaxies/"+tok[0];
      CVector pos( strtod(tok[1].c_str(),0), strtod(tok[2].c_str(),0), strtod(tok[3].c_str(),0) );
      CVector orientation( strtod(tok[4].c_str(),0), strtod(tok[5].c_str(),0), strtod(tok[6].c_str(),0) );
      double magnitude = strtod(tok[7].c_str(),0);
      orientation = orientation.Normalize();
      AddGalaxy(name, pos, orientation, magnitude);
    }
  }
}


void CRasterizer::setupCamera(int idx) {
  int i = idx%(int)RP.size;
  int j = (idx-i)/(int)RP.size;
  CVector p((((double)i-RP.size/2.0))/(double)RP.size,((double)j-RP.size/2.0)/(double)RP.size,0);
  CVector up = RP.up;
  if (up.Cross(RP.camera).Length()==0)
    up = CVector(0,0,1);
  
  CVector viewDir = (RP.target-RP.camera).Normalize();
  
  CVector left = viewDir.Cross(up).Normalize();
  up = left.Cross(viewDir).Normalize();
  up = up *p.x;
  left = left * p.y;
  
  RP.direction  = (viewDir+(up + left*RP.cameraAspectRatio)*RP.cameraWidth).Normalize();

}



void CRasterizer::setupLensPlane() {
  lensplane.Create(RP.size, RP.size);

  double S = 3.0;

  /*  for (int i=0;i<10;i++) {

    int i = rand()%(int)RP.size;
    int j = rand()%(int)RP.size;

    lensplane.putSmallCircle(i, j, RP.size/S, 10);
    }*/
  lensplane.putSmallCircle(RP.size/2, RP.size/2, RP.size/S, 512.0/RP.size);

  
  
}

void CRasterizer::setupLensing(int idx) {

  // Now perturb camera!

  int x = idx%(int)RP.size;
  int y = (idx-x)/(int)RP.size;
  

  // CVector lensPoint(float strength, float radius, int x, int y);
  CVector alpha = lensplane.lensPoint(1.0, 0.0, x, y);
  
  CVector xx = RP.up.Cross(RP.direction);
  CVector yy = xx.Cross(RP.direction);
  xx.Normalize();
  yy.Normalize();

  //  RP.direction =RP.direction + alpha*0.08;
  RP.direction =RP.direction + (xx*alpha.y + yy*alpha.x)*0.08;


}


void CRasterizer::renderPixel(int idx) {
      CVector isp1, isp2;

      double I = 0;
      for (int i=0;i<galaxies.size();i++) {
	
	CGalaxy* g = galaxies[i];
	
	bool is = CVector::intersectSphere(RP.camera-g->params.position, RP.direction, g->params.axis, isp1, isp2);
	float dist = (isp1-isp2).Length();

	
	if ((RP.camera-isp1).Length()>(RP.camera-isp2).Length()) {
	  CVector tmp = isp1;
	  isp1 = isp2;
	  isp2 = tmp;
	}

	CVector ispDisk;
	CVector::rayIntersectPlane(g->params.position, g->params.orientation, RP.camera, RP.direction, ispDisk);
	//  cout << ispDisk << endl;
	
	
	if (is) {
	  getIntensity(g, I, isp1, isp2);
	}
      }	

      image.Data[idx].height = I;
}

double CRasterizer::getIntensity(CGalaxy* g, double& I, CVector& isp1, CVector& isp2) {
  CVector origin = isp1;
  float length = (isp1-isp2).Length();
  CVector dir = (isp1-isp2).Normalize();
  
  int N = length/RP.rayStep;
  double scale = 1.0;
  float height = 0;
  
  
  CVector P;
  float val = 0;

  for (int i=0;i<N;i++) {
    CVector p = origin + dir*RP.rayStep*i;
    
    for (int j=0;j<g->components.size();j++) {
      I = g->components[j]->calculateIntensity(I, p, scale);
    }
    if (I<0)
      I = 0;
  }
}



double CRasterizer::getZfromFWHM(CVector galpos, double FWHM) {

  // first, set up two rays : one towards the centre, the other towards the viewport side

  FWHM /= (double)RP.size; // between zero and one


  CVector viewDir = (RP.target-RP.camera).Normalize();
  int i = 0;
  int j = RP.size/2;
  CVector p((((double)i-RP.size/2.0))/(double)RP.size,((double)j-RP.size/2.0)/(double)RP.size,0);
  CVector up = RP.up;
  if (up.Cross(viewDir).Length()==0)
    up = CVector(0,0,1);
  
  
  CVector left = viewDir.Cross(up).Normalize();
  up = left.Cross(viewDir).Normalize();



  up = up *p.x;
  left = left * p.y;


  CVector oRay = (viewDir.Normalize());
  CVector sRay  = (viewDir+(up + left*RP.cameraAspectRatio)*RP.cameraWidth).Normalize();


  double angle = sRay.Angle(oRay);
  
  // but what is x: 
  //  double fwhmscale = 0.45; // for example

  double fwhmscale = 0.26; // for example


  //  FWHM = 0.5;
  double x = 1 / (FWHM/fwhmscale);
  double z = x / tan(angle);

  return z;
}




double CRasterizer::renderAxis(CVector isp1, CVector isp2, CVector ispDisk) {
  double intensity = 0;
  float dist = (isp1-isp2).Length();
  double in = 1000;
  // draw curves and lines
  if (dist<0.2) {
    intensity = in;
  }
  double distDisk1 = (isp1-ispDisk).Length();
  double distDisk2 = (isp2-ispDisk).Length();
  if (distDisk1<0.01 || distDisk2 < 0.01) {
    intensity = in;
  }
  double d1 = CVector::pointToLineDistance(CVector(0,0,0), CVector(1,0,1).Normalize(), ispDisk);
  double d2 = CVector::pointToLineDistance(CVector(0,0,0), CVector(-1,0,1).Normalize(), ispDisk);
  if (d1<0.005 || d2<0.005) {
    intensity = in;
  }
  return intensity;
}




void CRasterizer::loadSpectra() {
  CGraph* g = new CGraph();

  try {
  double sm = 100;
  double sc = 1.0/1250.0;

  int N = 100;


  /*    CGraph stars, disk, bulge;
  int NN = 10000;
  stars.Initialize(NN);
  disk.Initialize(NN);
  bulge.Initialize(NN);
  for (int i=0;i<NN;i++) {
    stars.Index[i] = i;
    disk.Index[i] = i;
    bulge.Index[i] = i;

    float x = i/100.0;
    stars.Val[NN-i-1] = x*x/10;
    x-=35;
    if (x<0) x = 0;
    disk.Val[i] = x*x/10;
    disk.Val[i] = sqrt(x)*100;
    int j = NN - i - 1;
    bulge.Val[i] = 100  + 600*exp(-i/2500.0);////250 + j/25.0;

  }

  stars.Scale(1.0);
  disk.Scale(1.0);
  bulge.Scale(1.0);
  
  stars.SaveText("spectra/stars.txt");
  disk.SaveText("spectra/red.txt");
  bulge.SaveText("spectra/dust.txt");
  */
  
  //g->LoadText("spectra/blue.txt");
  g->LoadText("spectra/disk.lis");
  //g->Smooth(sm);
  g->Mul(sc);
  g->Scale(1.0);
  g->Crush(N);
  spectra.push_back(g);

  CGraph* g2 = new CGraph();

  g2->LoadText("spectra/red.txt");
  //g2->LoadText("spectra/Sbc.lis");
  //g2->Smooth(sm);
  g2->Mul(sc);
  g2->Scale(1.0);
  g2->Crush(N);
  spectra.push_back(g2);

  CGraph* g3 = new CGraph();

  //g3->LoadText("spectra/bulge.txt");
  g3->LoadText("spectra/bulge.lis");
  //g3->Smooth(sm);
  g3->Mul(sc);
  g3->Scale(1.0);
  g3->Crush(N);
  spectra.push_back(g3);


  CGraph* g4 = new CGraph();

  //  g3->LoadText("spectra/bulge.txt");
  g4->LoadText("spectra/dust.txt");
  g4->Smooth(sm);
  //g4->InvertX();


  g4->SaveText("spectra/dust2.txt");

  g4->Mul(sc);
  g4->Scale(1.0);
  g4->Crush(N);

  spectra.push_back(g4);



  }
  catch (string s) {
    cout << s << endl;
    exit(1);
  }

}
